System and method for tracking a global shape of an object in motion

ABSTRACT

A system and method for tracking a global shape of an object in motion is disclosed. One or more control points along an initial contour of the global shape are defined. Each of the one or more control points is tracked as the object is in motion. Uncertainty of a location of a control point in motion is represented using a number of techniques. The uncertainty to constrain the global shape is exploited using a prior shape model. In an alternative embodiment, multiple appearance models are built for each control point and the motion vectors produced by each model are combined in order to track the shape of the object. The movement of the shape of the object can be visually tracked using a display and color vectors.

CROSS REFERENCE TO RELATED APPLICATION

This application claims the benefit of U.S. Provisional Application Ser.No. 60/452,669, filed on Mar. 7, 2003, U.S. Provisional Application Ser.No. 60/473,425, filed on May 27, 2003, and U.S. Provisional ApplicationSer. No. 60/508,367 filed on Oct. 3, 2003 which are incorporated byreference in their entirety.

FIELD OF THE INVENTION

The present invention is directed to a method for tracking a shape thatis in motion, and more particularly, to a method for tracking a shapehaving linear constraints in the presence of heteroscedastic noise.

BACKGROUND OF THE INVENTION

For most visual tracking applications, measurement data are uncertainand sometimes missing: images are taken with noise and distortion, whileocclusions can render part of the object-of-interest unobservable.Uncertainty can be globally uniform; but in most real-world scenarios,it is heteroscedastic in nature, i.e., both anisotropic andinhomogeneous. A good example is the echocardiogram (ultrasound heartdata). Ultrasound is prone to reflection artifacts, e.g., specularreflectors, such as those that come from membranes. Because of thesingle “view direction”, the perpendicular surface of a specularstructure produces strong echoes, but tilted or “off-axis” surfaces mayproduce weak echoes, or no echoes at all (acoustic “drop out”). For anechocardiogram, the drop-out can occur at the area of the heart wherethe tissue surface is parallel to the ultrasound beam.

Due to its availability, relative low cost, and noninvasiveness, cardiacultrasound images are widely used for assessing cardiac functions. Inparticular, the analysis of ventricle motion is an efficient way toevaluate the degree of ischemia and infarction. Segmentation ordetection of the endocardium wall is the first step towardsquantification of elasticity and contractibility of the left ventricle.Examples of some existing methods include pixel-basedsegmentation/clustering approaches (e.g., Color Kinesis), variants ofoptical flow, deformable templates and Markov random process/fields, andactive contours/snakes. Some methods are employed in 2-Dimensional,3-Dimensional or 4-Dimensional (3D+time) space.

However, most existing segmentation or detection methods do not attemptto recover accurate regional motions of the endocardial wall, and inmost cases, motion components along the wall are ignored. Thissimplified treatment is also employed by contour trackers that searchonly along the normals of the current contour. This is not suitable forregional wall abnormality detection, because regional motion of anabnormal left ventricle is likely to be off the normal of the contour,not to mention that global motion, such as translation or rotation (dueto the sonographer's hand motion or respiratory motion the patient),causes off-normal local motion on the contour as well. It is desirableto track the global shape of endocardial wall as well as its localmotion, for the detection of regional wall motion abnormalities. Thisinformation can be used for further diagnosis of ischemia andinfarction.

In general, covariances can be assigned to image features or flowestimates that reflect underlying heteroscedastic noise. When the datais clean with a low overall noise level, the heteroscedastic nature maybe ignorable, and a global uncertainty can be substituted for the localestimates. However, for very noisy inputs, especially those withspatially varying structural noise, the information encoded in the localcovariance matrix becomes critical in ensuring reliable and robustinference of objects or underlying image structures.

It is a common practice to impose model constraints in a trackingframework. Examples include simple models such as blobs or parameterizedellipses, and complex models such as discriminative templates. In mostpractical cases, a subspace model is suitable for shape tracking, sincethe number of modes capturing the major shape variations is limited andusually much smaller than the original number of feature components usedto describe the shape. Furthermore, a Principal Component Analysis(PCA)-based eigenshape subspace can capture arbitrarily complicatedshape variations, which, in the original space, even with a very simpleparametric model, are highly nonlinear.

If a measurement vector is affected by heteroscedastic noise, anorthogonal projection into the constraining subspace is not onlyunjustified, but also very damaging in terms of information loss. It canonly be justified for the special case when the noise is both isotropicand homogeneous.

However, most existing work on subspace-constrained tracking does nottake into account the heteroscedastic noise in the measurements. In the“Point Distribution Model” or “Active Shape Model”, a PCA-based subspaceshape model is derived based on training shapes with landmark pointcorrespondence. The resulting subspace of eigenshapes captures the mostsignificant variations in the training data set. At detection time, amodel is perturbed to create synthetic images for matching against thetesting image at a candidate location. However, the measurement noisewas not modeled in this process.

Even when heteroscedastic noise characteristics are available, they weretypically disregarded during the subspace model fitting. For example, inone known approach where full covariance matrix was captured in themeasurements, a rather ad hoc thresholding is applied so that themeasurement mean is confined to a hyper-ellipsoid constraint defined bythe model covariance. This operation is independent of the measurementnoise.

Another known approach applies a two-step approach to impose a shapespace constraint in a Kalman filtering framework. The shape space is alinearly transformed affine subspace or eigen-subspace. However, theprojection into the shape space is orthogonal, without taking intoaccount the heteroscedastic noise of the measurement. Therefore, thisapproach leads to information loss during the projection.

Another known approach uses a Gaussian distribution to adaptively modelthe appearance of the object of interest (face in their case), which islearned using the EM algorithm. As in the present invention, localuncertainty is captured in the covariance matrix. The difference is thatthe present invention specifically studies the subspace modelconstraints and the critical choice of intersection over projection whenanisotropic uncertainty is present.

Another known approach uses a subspace constraint implicitly during theoptical flow estimation, which also utilizes flow uncertainties.Although in a different framework for a different application, presentinvention recognizes that “more reliable flow-vectors will have moreinfluence in the subspace projection process.”

Another known approach applies heteroscedastic regression for fittingellipses and fundamental matrices. The fitting is achieved in theoriginal space with parameterized models. In the present invention,parameterization of shape variations is avoided—it can be verycomplicated and highly nonlinear. Instead, the present invention buildssubspace linear probabilistic models through, e.g., PCA, and obtainclosed-form solutions on both the mean and covariance of the fitteddata.

Robust model matching relying on M-estimators or RANSAC has also beenapplied to limit or eliminate the influence of data components that areoutliers with respect to the model. Again, the locally (in space ortime) varying uncertainties are not exploited in these frameworks.

Other related approaches include data imputation, the practice of“filling in” missing data with plausible values. Work in this area isrooted in statistics with broad applications toward speech recognition,medical image analysis, and social science, etc. However, theformulation of data imputation problems typically assumes 0-1availability, i.e., a data component is either missing or available.There is a need for a unified framework for fusing subspace modelconstraints with information about the shape dynamic and theheteroscedastic nature of the measurement noise and about the shapedynamics.

SUMMARY OF THE INVENTION

The present invention is directed to a system and method for tracking aglobal shape of an object in motion. One or more control points along aninitial contour of the global shape are defined. Each of the one or morecontrol points is tracked as the object is in motion. Uncertainty of alocation of a control point in motion is estimated. One form torepresent the uncertainty is a covariance matrix. When employing asubspace shape constraint model, the uncertainty is exploited using anon-orthogonal projection and/or information fusion. Each subsequentcontour is displayed.

The present invention is also directed to a system for visually trackingmovement of a shape of an object. One or more first color vectors aregenerated to represent contraction of control points along the contourof the shape. One or more second control vectors are generated torepresent dilation of control points along the contour of the shape. Thefirst and second color vectors are periodically displayed therebydisplaying movement of the shape.

The present invention is also directed to a method for tracking a globalshape of an object in motion. One or more control points are definedalong the global shape. Each of the one or more control points aretracked as the object is in motion. Multiple appearance models are builtfor each control point. Motion vectors produced by each model fortracking the shape are combined.

BRIEF DESCRIPTION OF THE DRAWINGS

Preferred embodiments of the present invention will be described belowin more detail, wherein like reference numerals indicate like elements,with reference to the accompanying drawings:

FIG. 1 is a block diagram of an exemplary system for implementing amethod for shape tracking in accordance with the present invention;

FIG. 2 is an echocardiography image of a heart that illustrates areas ofacoustic drop out and estimated local wall motion uncertainties;

FIG. 3 is an echocardiography image of a left ventricle that illustratesan endocardial contour with localization uncertainties of its controlpoints that are strongly anisotropic and inhomogeneous;

FIG. 4 illustrates examples of an incremental PCA model and a stronglyadapted PCA (SA-PCA) model with different α values;

FIGS. 5 a-5 c illustrate echocardiography images of a left ventricle inwhich the endocardial wall is initialized (a) and tracked using IPCA (b)and SA-PCA (c);

FIGS. 6 a and 6 b illustrate echocardiography images of a left ventriclewherein the movement of the endocardial wall is tracked in accordancewith the present invention;

FIGS. 7 a-7 d illustrate test contours representing echocardiographyimages of an endocardial wall of a left ventricle from an apicalfour-chamber view in accordance with the present invention;

FIGS. 8 a-8 d illustrate test contours representing echocardiographyimages of an endocardial wall of a left ventricle from a parasternalshort axis view in accordance with the present invention;

FIG. 9 illustrates echocardiography images which exemplify a method forvisualizing the images in accordance with the present invention; and

FIG. 10 illustrates a multi-model control point based tracker inaccordance with an embodiment of the present invention.

DETAILED DESCRIPTION

The present invention is directed to a method for tracking shapes withlinear constraints in the presence of heterodastic noise. An examplewhere such a method would be utilized is for tracking the global shapeof a myocardial wall as well as its local motion to detect regional wallmotion abnormalities in the heart. The method may also be used to trackthe endocardial wall or epicardial wall of the heart. It is to beunderstood by those skilled in the art that the present invention may beused in other applications where shape tracking is useful such as, butnot limited to, recognizing movement of human features such as headmovements, facial features, hand movements or other body movements. Thepresent invention can also be used in 2 dimensional, 3 dimensional and 4dimensional (3D+time) medical analyses of anatomical structures such asthe heart, lungs or tumors that are evolving over time.

For purposes of describing the present invention, an example will bedescribed for tracking the endocardial wall of the left ventricle. FIG.1 illustrates an exemplary architecture of an echocardiograph systemthat uses a method for tracking the shape of an endocardial wall of aleft ventricle in accordance with the present invention. A medicalsensor 102, such as an ultrasound transducer is used to perform anexamination on a patient. The sensor 102 is used to obtain medicalmeasurements consistent with a particular medical examination. Forexample, a patient experiencing heart problems may have anechocardiogram performed to help diagnose the particular heart ailment.An ultrasound system provides two-, three-, and four(3D+time)-dimensional images of the heart from various perspectives.

The information obtained by the sensor 102 is communicated to aprocessor 104 which may be a workstation or personal computer. Theprocessor 104 converts the sensor data into an image that iscommunicated to display 108. The display 108 may also communicate othergraphical information or tables of information relating to the image. Inaccordance with the present invention, the processor 104 is alsoprovided with data representing an initial contour of the endocardialwall. The data may be provided manually by a user such as a physician orsonographer, or automatically by the processor 104. The contourcomprises a series of individual points, the movement of which istracked by the processor 104 and illustrated on display 108. Thespecifics regarding how the individual points are tracked will bedescribed in greater detail hereinafter.

In addition to data from the medical sensor 102, the processor 104 mayalso receive other data inputs. For example, the processor may receivedata from a database 106 associated with the processor 104. Such datamay include subspace models that represent potential contour shapes forthe endocardial wall. These subspace models may be images of leftventricles that are representative of a plurality of patients or may becomputer generated models of contour shapes based on statisticalinformation. The processor 104 tracks the individual points of thecontour shape using known approaches such as Bayesian kernel matching oroptical flow-based methods. Error accumulation during tracking isremedied by using a multi-template adaptive matching framework.Uncertainty of tracking is represented at each point in the form of acovariance matrix, which is subsequently fully exploited by a subspaceshape constraint using a non-orthogonal projection.

FIG. 2 illustrates a typical echocardiograph image of a heart The partof the endocardial wall of the left ventricle that has acoustic drop outis marked by solid ellipse 208. Estimations of local wall motion aredepicted by the doffed ellipses 202, 204 and 206. Because of theacoustic drop out, the endocardium wall is not always at the strongestedge in the image. In accordance with the present invention, a unifiedframework is used for fusing subspace model constraints with informationabout the shape dynamic and the heteroscedastic nature of themeasurement noise and about the shape dynamics. The subspace model cantake the form of a specific subspace distribution, e.g., a Gaussian, ora simple subspace constraint, e.g., the eigenspace model.

The present invention tracks individual control points on a contour thatrepresents the endocardium wall. The tracking may be performed using aBayesian kernel matching approach or a flow-based approach. An exampleof a Bayesian kernel matching approach is described in the articleauthored by co-inventor, Dorin Comaniciu and titled: Bayesian KernelTracking, Annual Conf. of the German Society for Pattern Recognition(DAGM'02), Zurich, Switzerland, 438-445, 2002 which is incorporated byreference in its entirety. An example of an optical flow based methodfor tracking the individual points is described in co-pendingapplication Ser. No. 10/681,702 entitled “Density Estimation-BasedInformation Fusion for Multiple Motion Computation” which isincorporated by reference in its entirety. In accordance with thepresent invention, error accumulation during tracking is remedied usinga multi-template adaptive matching framework, taking advantage of theperiodic nature of cardiac motion. Uncertainty of tracking isrepresented at each point in the form of a covariance matrix, which issubsequently fully exploited by a subspace shape constraint using anon-orthogonal projection.

The tracking framework is a two-step, iterative process over the imagesequences. An initial contour with control points is drawn on the firstimage (either automatically or manually); then, for every subsequentimage in order, every control point is first tracked independently, withanisotropic uncertainty also recorded. As the second step, the newcontour is projected into a feasible subspace using non-orthogonalprojection. The feasible subspace is learned based training contours andalso adapted to the current case using the initial contour available tothe tracker, also taking into account the confidence of initialization(i.e., a manual initialization carries a higher confidence than a fullyautomatic one).

As indicated above, during the tracking process multiple templates areemployed. The use of the multiple templates results in a more accuraterepresentation of the shape statistics for the current case. In theBayesian kernel matching approach, a first template is taken from theinitialized frame. Subsequent templates are added when they are bothsufficiently different from the existing templates, and alsosufficiently informative for localization purpose, which is measured bykernel matching with itself: a more informative patch is one that has ahigher confidence when matched to itself.

The decision to use more than one templates is based on the observationthat the heart motion is periodic, therefore, different appearancepatterns within one cycle will all reappear in later cycles. To matchwith multiple templates, match is made to each template and the one withthe best match is selected. Or, to save computation, only theneighboring templates of the previous matching template may be matched,exploiting again the periodic motion.

Uncertainty on matching location is calculated in the neighborhood ofthe optimal location, with the likelihood map calculated in the same wayusing the kernel matching process. The likelihood surface is then usedfor estimating a covariance matrix, through, e.g., fitting with atwo-dimensional Gaussian distribution or the inversion of a weightedestimation of the Hessian matrix.

In the case of the flow-based approach, a typical optical flowimplementation uses only neighboring frames, making the tracking processover long sequences susceptible to error accumulation and drifting.Multiple templates are employed for flow calculation, much the same wayas for kernel matching. A new template is added whenever the flowuncertainty is high, while the local gradient is distinctive fromexisting templates.

After the new locations are obtained for every control point, the nextstep is to constrain the overall shape by a statistical model thatcaptures the “legitimate” shape variations of a human heart. PCA-basedshape models, which are also known as point distribution models, oractive shape models are used.

Because the uncertainties over the heart contour are neitherhomogeneous, some areas are worse than others due to, say, signaldropouts; nor isotropic, for example, there is worse localization alongan edge than along the gradient direction. FIG. 3 shows an example ofsuch anisotropic and inhomogeneous noise across a contour. As can beseen, individual points, such as points 302, 304, 306 and 308 have beeninitially identified along the contour of the endocardial wall. Acertainty measurement is taken for each point and an ellipsoid iscreated around each point which indicates the level of certainty thatthat particular point is in the correct location; the larger theellipsoid, the greater the level of uncertainty as to the location of aparticular point. As seen in FIG. 3, point 302 has a relatively smallellipsoid surrounding it, thereby indicating that there is a high degreeof certainty as to the location of the point on the contour. Points 306and 308 have relatively large ellipsoids surrounding them, therebyindicating a large degree of uncertainty as to the location of thesepoints.

Taking into account a general uncertainty model over the new contourpoints, the optimal projection will no longer be orthogonal. The optimalsolution in the subspace is in fact the maximal likelihood shape on theintersected distribution in the shape model subspace. Furthermore, ifthere is a distribution model in the subspace, there is no reason forignoring such extra information. In the following the detailed analysisis provided as to how to perform non-orthogonal projection; and, if thesubspace model distribution is available, how to fuse the information ofsuch a model with the uncertain input.

The present invention uses a subspace model-based fusion approach. Giventwo noisy measurements of the same n-dimensional variable x, eachcharacterized by a multidimensional Gaussian distribution, p₁ and p₂,the maximum likelihood estimate of x is the point with the minimal sumof Mahalanobis distances to the two centroids: x*=argmin d_(m) withd _(m)(x−x ₁)^(T) C ₁ ⁻¹(x−x ₁)+(x−x ₂)^(T) C ₂ ⁻¹(x−x ²)  (1)

Taking derivative with respect to x, we get:x*=C(C ₁ ⁻¹ x+C ₂ ⁻¹ x ₂)C=(C ₁ ⁻¹ +C ₂ ⁻¹)⁻¹  (2)

This is also known as the best linear unbiased estimate (BLUE) of x.

Assume that one of the Gaussians is in a subspace of dimension p, e.g.,C₂ is singular. With the singular value decomposition of C₂=UΛU^(T),where U=[u₁, u₂, . . . , u_(n)], with u_(i)'s orthonormal and Λ=diag{λ₁,λ₂, . . . , λ_(p), 0, . . . , 0}, we rewrite Mahalanobis distance to x₂in Eq. (1) in the canonical form:

$\begin{matrix}{d_{m,2} = {{\left( {x - x_{2}} \right)^{T}{C_{2}^{- 1}\left( {x - x_{2}} \right)}} = {\sum\limits_{i = 1}^{n}{\lambda_{i}^{- 1}\left\lbrack {U^{T}\left( {x - x_{2}} \right)} \right\rbrack}^{2}}}} & (3)\end{matrix}$

When λ_(i) tends to 0, d_(m,2) goes to infinity, unless U₀ ^(T)x=0,where U₀=[u_(p+1), u_(p+2), . . . , u_(n)]. (Here we have assumed,without loss of generality, the subspace passes through the origin ofthe original space. Since x₂ resides in the subspace, U₀ ^(T)x₂=0.)

Therefore, we only need to retain in Eq. (3) those terms correspondingto a non-zero λ_(i):

$\begin{matrix}{{d_{m,2} = {{\sum\limits_{i = 1}^{p}{\lambda_{i}^{- 1}\left\lbrack {U_{p}^{T}\left( {x - x_{2}} \right)} \right\rbrack}^{2}} \equiv {\left( {x - x_{2}} \right)^{T}{C_{2}^{+}\left( {x - x_{2}} \right)}}}},} & (4)\end{matrix}$

where C⁺ is the pseudoinverse of C, and U_(p)=[u₁, u₂, . . . , u_(p)].

Furthermore, because U₀ ^(T)x=0, x can be expressed in another form toreflect this constraint:x=UU ^(T) x=U([U _(p) |U ₀]^(T) x)=U[y|1]=U _(p) y,  (5)

for a 1×p vector y. Eq. (1) now takes the following general form despitesingularity:d _(m)=(U _(p) y−x ₁)^(T) C ₁ ⁻¹(U _(p) y−x ₁)+(U _(p) y−x ₂)^(T) C ₂⁺(U _(p) y−x ₂)  (6)

Taking derivative with respect to y yieldsy*=C _(y*) U _(p) ^(T)(C ₁ ⁻¹ x ₁ +C ₂ ⁺ x ₂)  (7)C _(y*) =[U _(p) ^(T)(c ₁ ⁻¹ +C ₂ ⁺)U _(p)]⁻  (8)x*=U _(p) y*=C _(x)*(C ₁ ⁻¹ x ₁ +C ₂ ⁺ x ₂)  (9)C _(x*) =U _(p) C _(y*) U _(p) ^(T)  (10)

Eq. (7) indicates a fusion in the information space, followed by atransformation into the subspace. Eq. (9) is a coordinate transform backinto the original space. It can be shown that C_(x*) and C_(y*) are thecorresponding covariance matrices for x* and y*.

Notice that this solution is not a simple generalization of Eq. (2) bysubstituting pseudoinverses for regular inverses, which will notconstrain x* to be in the subspace.

Alternatively, we can write out Eq. (7) and (8) asy*(U _(p) ^(T) C ₁ ⁻¹ U _(p)+Λ_(p) ⁻¹)⁻¹(U _(p) ^(T) C ₁ ⁻¹ x ₁ +Λ _(p)⁻¹ y ₂)  (11)

Here y₂ is the transformed coordinates of x₂ in the subspace spanned byU_(p). and Λ_(p)=diag{λ₁, λ₂, . . . , λ_(p)}.

Eq. (11) can also be shown as a BLUE fusion in the subspace of twodistributions, one is N(y₂, Λ_(p)) and the other is the “intersection”(not projection!) of N(x₁, C₁) in the subspace, N((U_(p) ^(T)C₁⁻¹U_(p))⁻¹U_(p) ^(T)C₁ ⁻¹x₁, (U_(p) ^(T)C₁ ⁻¹U_(p))⁻¹).

The use of a statistical shape model learned from a large pool oftraining samples to guide the contours from a specific heart cansometimes be problematic. In accordance with the present invention, thetraining samples are used to obtain a shape model of the current heart,and not a generic heart. Therefore, there is a strong motivation toadapt the generic model toward what is known for the current case. Inaccordance with the present invention, the initial contour (manual orthrough automatic detection) of the endocardial wall of the heart of thepatient is examined to adapt the existing PCA model.

In determining the actual contour of the endocardial wall and trackingits movement, two approaches are considered: one with the initialcontour being deterministic, and the other with the initial contourbeing uncertain (this can be the case when the initial contour is froman automatic detection algorithm, which also provides uncertainties).

When the initial contour is assumed to be a point (deterministicallycertain), a strongly-adapted-PCA (SA-PCA) model is used to track themovement of the point. It is assumed that the old PCA model (excludingthe current case) and the initialized contour for the current casejointly represent the variations of the current case, but with relativeenergy (i.e., representative power) being α and (1−α), respectively,with 0<α<1. In other words, we assume that a portion of shape variationsof the current case will be represented in the generic model, while therest is captured in the direction of the initial contour.

The PCA model is denoted by its mean, eigenvalue matrix, and eigenvectormatrix by x_(m), Λ, and U, respectively. If the original full covariancematrix C is stored (this would be the case when the originaldimensionality is not forbiddingly high), the adapted mean andcovariance matrix are simply the weighted sum of the two contributingsources:x _(m,new) =αx _(m)+(1−α)_(ix)  (12)

$\begin{matrix}\begin{matrix}{C_{new} = {{\alpha\left( {C + {\left( {x_{m} - x_{m,{new}}} \right)\left( {x_{m} - x_{m,{new}}} \right)^{T}}} \right)} +}} \\{\left( {1 - \alpha} \right)\left( {x - x_{m,{new}}} \right)\left( {x - x_{m,{new}}} \right)^{T}} \\{= {{\alpha\; C} + {{\alpha\left( {1 - \alpha} \right)}\left( {x - x_{m}} \right)\left( {x - x_{m}} \right)^{T}}}}\end{matrix} & (13)\end{matrix}$

Eigenanalysis can be performed on C_(new) to obtain the new subspacemodel.

In the case that C is not stored but only {x_(m), Λ, U} is available inthe subspace, through straight algebraic manipulations we can arrive atthe adapted eigenanalysis results {x_(m,new), Λ_(new), U_(new)} asfollows:

The initial contour x has a subspace component as x_(s)=U^(T)x_(d),where x_(d)=x−x_(m), and a residual vector as x_(r)=(x−x_(m))−U x_(s).Let x_(ru) be the normalized version of x_(r) with norm 1 (or zero ifx_(r) is zero norm).

The adapted eigenvector matrix that represents the combined energy willhave the following form:U _(new) =[U, x _(ru) ]R  (14)

R and Λ_(new) will be the solutions to the following eigenanalysisproblem:

$\begin{matrix}{{\left( {{\alpha\begin{bmatrix}\Lambda & 0 \\0^{T} & 0\end{bmatrix}} + {{\alpha\left( {1 - \alpha} \right)}\begin{bmatrix}{x_{s}x_{s}^{T}} & {e_{r}x_{s}} \\{e_{r}x_{s}^{T}} & e_{r}^{2}\end{bmatrix}}} \right)R} = {R\;\Lambda_{new}}} & (15)\end{matrix}$

where e_(r)=x_(ru) ^(T)(x−x_(m)) is the residual energy.

Above formulae are more general than IPCA, with tunable energy ratiosbetween the new data and the old data. These become equivalent to IPCAif we set α to be the fraction of points in the model versus the totalnumber of points. Typically, this will be a number very close to 1 sincethe number of contours in the training set is usually large. With a setat a smaller value (say, 0.5), the PCA model is strongly adapted towardthe current case, hence the name. FIG. 4 shows a simple two-dimensionalillustration of IPCA and SA-PCA with different α values.

Point 404 represents the current case. Each x 402 represents a trainingpoint corresponding to a particular model. Ellipse 406 shows theoriginal model distribution. Ellipse 408 illustrates an Incremental PCAmodel which corresponds to a strongly adapted PCA model in which α=0.99.Ellipse 410 illustrates a strongly adapted PCA model in which α=0.5.Ellipse 412 illustrates a strongly adapted PCA model in which α=0.1.Each ellipse depicts a 90% equal-probable contour of the correspondingdistribution.

Indeed, the contours from the current heart are much more likely toresemble the initial contour of the same heart than those contours inthe generic training set, especially if the current heart has anirregular shape that is not well represented in the training set. In oursystem, we set α at 0.5. This permits a strong influence from theinitial contour: 50% of model energy is from the initial contour (thisreduces to only 1% if IPCA is applied with a model trained with 99examples).

FIGS. 5 a-5 c show a comparison of IPCA and SA-PCA. This parasternalshort axis view has an irregular shape (with a concave part), while thetraining set is overwhelmingly populated with circular shapes. Theincremental PCA model, taking in the initial contour (FIG. 5 a) but witha very small weight (<0.01%), fails to capture the concave nature of thecurrent shape; and has constrained the contours to a typical circularshape (FIG. 5 b). This result is in fact the same as that obtained usingthe old PCA model without an incremental step. Our adaptive PCA modelwith α=0.5, fits much better to the true borders (FIG. 5 c).

A subtle while important interaction between the non-orthogonalprojection with fusion approach and the SA-PCA model adaptation is asfollows: the fusion with the model mean and covariance is necessary tofilter out contours that are in the subspace but too far from the modeldistribution. However, this stronger constraint (than orthogonalprojection) will inevitably alter low likelihood, or outlier heartcontours. The SA-PCA model mediates this dilemma through a strong modelshift toward the direction of the current case using the extrainformation provided in the initial contour given.

A Kalman filter fuses information from the prediction defined by adynamic process and from noisy measurements. When applied for shapetracking, additional global constraints are necessary to stabilize theoverall shape in a feasible range. The present invention utilizes aunified fusion framework to incorporate subspace model constraints intoa Kalman filter.

For a Kalman filter, the measurement update equation has the followingform:x_(k+1|k+1) =x _(k+1|k) +K(z _(k+1) −Hx _(k+1|k))  (16)whereK=P _(k+1|k) H ^(T)(HP _(k+1|k) H ^(T) +R)⁻¹  (17)P _(k+1|k+1)=(I−KH)P _(k+1|k)  (18)P _(k+1|k) =SP _(k|k) S ^(T) +Q  (19)

Here P is the estimation error covariance, and x_(i|j) is the stateestimate at time i given the state at time j. The measurement model isz_(k)=H x_(k)+r_(k), where r_(k) represents measurement noise withcovariance R. The system/process model is x_(k+1)=S x_(k)+q_(k), whereq_(k) represents system noise with covariance Q.

Using the approach above, the update equations of the Kalman filter withsubspace constraint and heteroscedastic noise are given byx _(k+1|k+1) =P _(k+1|k+1)((P _(k|k) +Q)⁺ x _(k|k) +R ⁻¹ z+C ₂ ⁺ x₂)  (20)P _(k+1|k+1) =U _(p) [U _(p) ^(T)(((P _(k|k) +Q)⁺ +R ⁻¹ +C ₂ ⁺)U _(p) ]U_(p) ^(T)  (21)where it is assumed that the system noise covariance Q is contained inthe subspace. Observe the symmetry of the solution which combines allthe available knowledge in the information space. These equationsprovide a unified fusion of the system dynamics, subspace constraint,and noise information. They represent the complete representation ofvarious uncertainties that affect the tracking system.

In accordance with an alternative embodiment of the present invention, atracking technique is used that relies on control point based objectrepresentations and on robust fusion to integrate model informationacross frames. A set of control point based object representations aremaintained and acquired at different time instances. The estimatedmotion suggested by the control points are robustly combined todetermine the next position of the object. Visual tracking of the objectis achieved by maintaining several models over time. The result is anonparametric representation of the probability density function thatcharacterizes the object appearance. Tracking is performed by obtainingindependently from each model a motion estimate and its uncertaintythrough optical flow. The final estimate for each control point iscomputed using a robust fusion technique such as a Variable-BandwidthDensity-based Fusion (VBDF) procedure. VBDF computes the location of themost significant mode of the displacements density function while takinginto account their uncertainty. The VBDF procedure manages the multipledata sources and outliers in the displacement estimates. Occlusions arenaturally handled through the estimate uncertainty for large residualerrors. The model is divided in several regions for which the flow iscomputed independently. The residual alignment error is used to computethe scale of the covariance matrix of the estimate, therefore reducingthe influence of the unreliable displacements.

Object tracking challenges due to occlusions and appearance variationsare handled in the present invention through a multi-modal control pointbased approach. Maintaining several representatives for a 2 dimensionalappearance model does not restrict it to a unimodal distribution and theVDBF mechanism robustly integrates multiple estimates to determine themost dominant motion for each control point. To model changes duringtracking, several exemplars of the object appearance over time aremaintained.

FIG. 10 illustrates an exemplary multi-model control point based trackerin accordance with the present invention. The top three frames 1002,1004 and 1006 illustrate the current exemplars in the model set, eachhaving associated a set of overlapping control points. A control pointbased approach is more robust than a global representation and thereforeless sensitive to illumination changes and pose. Another advantage ofthe present invention is that partial occlusion can be handled at thecontrol point level by analyzing the matching likelihood.

Each control point is processed independently; its location andcovariance matrix is estimated in the current image with respect to allof the model templates. For example, one of the control points isillustrated by gray rectangle 1014 and its location and uncertainty withrespect to each model is shown in I_(new) 1008. The VBDF robust fusionprocedure is applied to determine the most dominant motion (mode) withthe associated uncertainty as illustrated in frame 1010. Note thevariance in the estimated location of each control point due toocclusion or appearance change.

The location of the control points in the current frame is furtherconstrained by a global parametric motion model. A similaritytransformation model and its parameters are estimated using theconfidence in each control point location. Therefore, the reliablecontrol points contribute more to the global motion estimation. Thecurrent frame is added to the model set if the residual error so thereference appearances is relatively low. The threshold is chosen suchthat the images are not added where the object has significantocclusion. The number of templates in the model can be variable orfixed. If the number of templates are fixed, a scheme can be provided inwhich to discard certain templates (e.g., discard the oldest template).

A VBDF estimator is based on nonparametric density estimation withadaptive kernel bandwidths. The VBDF estimator is defined as thelocation of the most significant mode of the density function. The modecomputation is based on the variable-bandwidth mean shift technique in amultiscale optimization framework.

Let x_(i)εR^(d), i=1 . . . n be the available d-dimensional estimates,each having an associated uncertainty given by the covariance matrixC_(i). The most significant mode of the density function is determinediteratively in a multiscale fashion. A bandwidth matrix H_(i)=C_(i)+α²Iis associated with each point x_(i), where I is the identity matrix andthe parameter a determines the scale of the analysis. The sample pointdensity estimator at location x is defined by

$\begin{matrix}{{\hat{f}(x)} = {\frac{1}{{n\left( {2\;\pi} \right)}^{d/2}}{\sum\limits_{i = 1}^{n}{\exp\left( {{- \frac{1}{2}}{D^{2}\left( {x,x_{i},H_{i}} \right)}} \right)}}}} & (22)\end{matrix}$where D represents the Mahalanobis distance between x and x_(i)D ²(x,x _(i) , H _(i))=(x−x _(i))^(⊥) H _(i) ⁻¹(x−x _(i))  (23)

The variable bandwidth mean shift vector at location x is given by

$\begin{matrix}{{m(x)} = {{{H_{h}(x)}{\sum\limits_{i = 1}^{n}{{\omega_{i}(x)}H_{i}^{- 1}x_{i}}}} - x}} & (24)\end{matrix}$where H_(h) represents the harmonic mean of the bandwidth matricesweighted by the data-dependent weights ω_(i)(x)

$\begin{matrix}{{H_{h}(x)} = {\left( {\sum\limits_{i = 1}^{n}{{\omega_{i}(x)}H_{i}^{- 1}}} \right)^{- 1}.}} & (25)\end{matrix}$

The data dependent weights computed at the current location x have theexpression

$\begin{matrix}{{\omega_{i}(x)} = \frac{\frac{1}{{H_{i}}^{1/2}}{\exp\left( {{- \frac{1}{2}}{D^{2}\left( {x,x_{i},H_{i}} \right)}} \right)}}{\sum\limits_{i = 1}^{n}{\frac{1}{{H_{i}}^{1/2}}{\exp\left( {{- \frac{1}{2}}{D^{2}\left( {x,x_{i},H_{i}} \right)}} \right)}}}} & (26)\end{matrix}$

and note that they satisfy

${\sum\limits_{i = 1}^{n}\;{\omega_{i}(x)}} = 1.$

It can be shown that the density corresponding to the point x+m(x) isalways higher or equal to the one corresponding to x. Thereforeiteratively updating the current location using the mean shift vectoryields a hill-climbing procedure which converges to a stationary pointof the underlying density.

The VBDF estimator finds the most important mode by iteratively applyingthe adaptive mean shift procedure at several scales. It starts from alarge scale by choosing the parameter a large with respect to the spreadof the points x_(i). In this case the density surface is unimodaltherefore the determined mode will correspond to the globally densestregion. The procedure is repeated while reducing the value of theparameter a and starting the mean shift iterations from the modedetermined at the previous scale. For the final step the bandwidthmatrix associated to each point is equal to the covariance matrix, i.e.H_(i)=C_(i).

The VBDF estimator is a powerful tool for information fusion with theability to deal with multiple source models. This is important formotion estimation as points in a local neighborhood may exhibit multiplemotions. The most significant mode corresponds to the most relevantmotion.

Consider that we have n models M₀, M₁ . . . M_(n). For each image wemaintain c components with their location denoted by x_(ij), i=1 . . .n, j=1 . . . c. When a new image is available, the location and theuncertainty for each component and for each model is estimated. Thisstep can be done using several techniques such as ones based on imagecorrelation, spatial gradient or regularization of spatio-temporalenergy. Using the VBDF technique, the result is the motion estimatex_(ij) for each component and its uncertainty C_(ij). Thus x_(ij)represents the location estimate of component j with respect to model i.The scale of the covariance matrix is also estimated from the matchingresidual errors. This will increase the size of the covariance matrixwhen the respective control point is occluded therefore occlusions arehandled at the control point level.

The VBDF robust fusion technique is applied to determine the mostrelevant location x_(j) for component j in the current frame. The modetracking across scales results in

$\begin{matrix}\begin{matrix}{{\hat{x}}_{j} = {{C\left( {\hat{x}}_{j} \right)}{\sum\limits_{i = 1}^{n}{{\omega_{i}\left( {\hat{x}}_{j} \right)}{\hat{C}}_{ij}^{- 1}{\hat{x}}_{ij}}}}} \\{{C\left( {\hat{x}}_{j} \right)} = {\left( {\sum\limits_{i = 1}^{n}{{\omega_{i}\left( {\hat{x}}_{j} \right)}{\hat{C}}_{ij}^{- 1}}} \right)^{- 1}.}}\end{matrix} & (27)\end{matrix}$with the weights ω_(i) defined as in (26)

$\begin{matrix}{{\omega_{i}\left( {\hat{x}}_{j} \right)} = \frac{\frac{1}{{C_{ij}}^{1/2}}{\exp\left( {{- \frac{1}{2}}{D^{2}\left( {{\hat{x}}_{j},{\hat{x}}_{ij},C_{ij}} \right)}} \right)}}{\sum\limits_{i = 1}^{n}{\frac{1}{{C_{ij}}^{1/2}}{\exp\left( {{- \frac{1}{2}}{D^{2}\left( {{\hat{x}}_{j},{\hat{x}}_{ij},C_{ij}} \right)}} \right)}}}} & (28)\end{matrix}$and note that they satisfy

${\sum\limits_{i = 1}^{n}\;{\omega_{i}(x)}} = 1.$

It can be shown that the density corresponding to the point x+m(x) isalways higher or equal to the one corresponding to x. Thereforeiteratively updating the current location using the mean shift vectoryields a hill-climbing procedure which converges to a stationary pointof the underlying density.

The VBDF estimator finds the most important mode by iteratively applyingthe adaptive mean shift procedure at several scales. It starts from alarge scale by choosing the parameter a large with respect to the spreadof the points x_(i). In this case the density surface is unimodaltherefore the determined mode will correspond to the globally densestregion. The procedure is repeated while reducing the value of theparameter a and starting the mean shift iterations from the modedetermined at the previous scale. For the final step, the bandwidthmatrix associated to each point is equal to the covariance matrix, i.e.,H_(i)=C_(i).

The VBDF estimator is a powerful tool for information fusion with theability to deal with multiple source models. This is important formotion estimation as points in a local neighborhood may exhibit multiplemotions. The most significant mode corresponds to the most relevantmotion.

An example for tracking multiple control point models will now bedescribed in accordance with the present invention. There are n modelsM₀, M₁, . . . , M_(n). For each image c control points with theirlocation denoted by x_(ij), i=1 . . . n, j=1 . . . c are maintained.When a new image is available, the location and uncertainty for eachcontrol point and for each model are estimated. This step can be doneusing several techniques such as ones based on image correlation,spatial gradient or regularization of spatio-temporal energy. Using theVBDF technique, the result is the motion estimate x_(ij) for eachcontrol point and its uncertainty C_(ij). Thus x_(ij) represents thelocation estimate of component j with respect to model i. The scale ofthe covariance matrix is also estimated from the matching residualerrors. This will increase the size of the covariance matrix when therespective control point is occluded therefore occlusions are handled atthe control point level.

The VBDF robust fusion technique is applied to determine the mostrelevant location x_(j) for component j in the current frame. The modetracking across scales results in

$\begin{matrix}{{{\hat{x}}_{j} = {{C\left( {\hat{x}}_{j} \right)}{\sum\limits_{i = 1}^{n}{{\omega_{i}\left( {\hat{x}}_{j} \right)}{\hat{C}}_{ij}^{- 1}{\hat{x}}_{ij}}}}}{{C\left( {\hat{x}}_{j} \right)} = {\left( {\sum\limits_{i = 1}^{n}{{\omega_{i}\left( {\hat{x}}_{j} \right)}{\hat{C}}_{ij}^{- 1}}} \right)^{- 1}.}}} & (29)\end{matrix}$with the weights ω_(i) as defined in (28)

$\begin{matrix}{{\omega_{i}\left( {\hat{x}}_{j} \right)} = \frac{\frac{1}{{C_{ij}}^{1/2}}{\exp\left( {{- \frac{1}{2}}{D^{2}\left( {{\hat{x}}_{j},{\hat{x}}_{ij},C_{ij}} \right)}} \right)}}{\sum\limits_{i = 1}^{n}{\frac{1}{{C_{ij}}^{1/2}}{\exp\left( {{- \frac{1}{2}}{D^{2}\left( {{\hat{x}}_{j},{\hat{x}}_{ij},C_{ij}} \right)}} \right)}}}} & (30)\end{matrix}$

Following the computation of each control point, a weighted rectanglefitting is carried out with the weights given by the covariance matrixof the estimates. The image patches are related by a similaritytransform T defined by 4 parameters. The similarity transform of thedynamic control point location x is characterized by the followingequations.

$\begin{matrix}{{T(x)} = {{\begin{pmatrix}a & {- b} \\b & a\end{pmatrix}x} + \begin{pmatrix}t_{x} \\t_{y}\end{pmatrix}}} & (31)\end{matrix}$where t_(x), t_(y) are the translational parameters and a, bparameterize the 2D rotation and scaling.

The minimized criterion is the sum of Mahalanobis distances between thereference location x⁰ _(j) and the estimated ones x_(j) (j^(th) controlpoint location in the current frame).({circumflex over (x)} _(j) −T(x _(j) ⁰))^(T) C({circumflex over (x)}_(j))⁻¹({circumflex over (x)} _(j) −T(x _(j) ⁰)).  (32)

Minimization is done through standard weighted least squares. Note thatbecause the covariance matrix is used for each control point, theinfluence of points with high uncertainty is reduced.

After the rectangle is fitted to the tracked control points, the dynamiccomponent candidate inside the rectangle is uniformly resampled. It isassumed that the relative position of each control point with respect tothe rectangle does not change a lot. If the distance of the resampleposition and the track position computed by the optical flow of acertain control point is larger than a tolerable threshold, the trackposition is regarded as an outlier and replaced with a resampled point.The current image is added to the model set if sufficient control pointshave low residual error. The median residual error between the modelsand the current frame is compared with a predetermined threshold T_(h).

Given a set of models M₀, M₁, . . . , M_(n) in which the component j haslocation x_(ij) in frame i, our object tracking algorithm can besummarized by the following steps:

-   -   1. Given a new image I_(f) compute {circumflex over (x)}_(ij)        ^((f)) through robust optical flow starting from {circumflex        over (x)}_(j) ^((f−1)), the location estimated in the previous        frame;    -   2. For j=1 . . . c, estimate the location {circumflex over        (x)}_(j) ^((f)) of component j using VBDF estimator resulting in        Eq. (6);    -   3. Constrain the component location using the transform computed        by minimizing (9);    -   4. Add the new appearance to the model set if its median        residual error is less that T_(h).

The proposed multi-template framework can be directly applied in thecontext of shape tracking. If the tracked points represent the controlpoints of a shape modeled by splines, the use of the robust fusion ofmultiple position estimates increases the reliability of the locationestimate of the shape. It also results in smaller corrections when theshape space is limited by learned subspace constraints. If the contouris available, the model templates used for tracking can be selectedonline from the model set based on the distance between shapes.

FIGS. 6 a and 6 b illustrate echocardiography images of a left ventriclein which movement of the endocardial wall is tracked in accordance withthe present invention. Frame 1 for each FIG. 602, 610 illustrates theinitial contour of the wall as shown by the dots in the image. Thesubsequent frames 604, 606, 608, 612, 614, 616 illustrate how themovement of the wall is tracked over time. Measurements for each dot aremade in accordance with the present invention.

EXAMPLE

Hand traced heart contours in echocardiography images are used as thetraining set in this example. Both apical four-chamber (open contour)and parasternal short axis (closed contour) views are tested asillustrated in FIGS. 7 a-7 d and FIGS. 8 a-8 d. Landmark points are alsoannotated for each contour. FIG. 7 b and FIG. 8 b show the set oftraining contours drawn together. The eigenanalysis approach is appliedwhere every contour is a vector with coordinates of the ordered landmarkpoints as its components (34 components for the open contour and 36 forthe closed contour). Then PCA is performed on the matrix whose columnsare the training vectors. And the eigenvalues are used to form adiagonal matrix as the model covariance matrix in the subspace. Themeasurement covariance matrix is adjusted to test different scenarios.

FIGS. 7 a and 8 a illustrate test contours where the solid curves 704,802 are the ground truth and the dashed curves 702, 804 are the noisymeasurements. It is desired to find a contour in the aigen-contoursubspace that is the closest to the current noisy contour, with moreuncertain point adjusted more. FIGS. 7 b and 8 b depict the trainingsets 706 and 806 of contours. FIG. 7 c and FIG. 8 c show the result whenorthogonal projection is applied to the contours of FIG. 7 a and FIG. 8a (with isotropic covariance depicted by small circles 712, 812 aroundthe landmark points). Solid curves 710 and 808 are the ground truth anddashed curves 708 and 810 are the noisy measurements. FIG. 7 d and FIG.8 d show the results. Solid curves 714 and 814 are the ground truth anddashed curves 718 and 816 are the noisy measurements. Ellipses 716, 818and 820 are the equal probably contours depicting the covariance at eachcontrol point. It is seen that the result is much closer to the groundtruth. It is not a perfect fit because the training data is ratherdifferent from the test data, therefore, some small shape deformation inthe test contour may not be realizable in the subspace.

In accordance with the present invention, the contours can be viewedusing a visualization tool to facilitate the diagnosis. To make thedoctor diagnose the heart diseases from the echocardiograms more easily,LV contours are displayed by color vectors in real time. The length ofthe vectors represents the displacement magnitude while the direction ofthe vectors represents the movement direction of the points on the LVwall contours. For exemplary purposes, the color orange is used todepict when the motion is contraction and the color dark blue is used todepict when the motion is dilation. FIG. 9 shows one of the representingresults. The movements of the points on the LV contours are tracked. Themovements are smoothed by Gaussian in both temporal and spacious domain.To make it easier seen by the doctors, the length of the vectors ismagnified by three times.

By this visualization method, the doctor can easily see the motion ofeach segment of the endocardium with the magnitude and the directions.Combination of the global motion compensation with this visualizationmethod can make the doctors easily see the contraction direction andmagnitude in all the segments of the LV. An ischemic region or otherabnormal regions can be identified directly and easily by human eyes.

Having described embodiments for method for determining featuresensitivity during a medical examination, it is noted that modificationsand variations can be made by persons skilled in the art in light of theabove teachings. It is therefore to be understood that changes may bemade in the particular embodiments of the invention disclosed which arewithin the scope and spirit of the invention as defined by the appendedclaims. Having thus described the invention with the details andparticularity required by the patent laws, what is claimed and desiredprotected by Letters Patent is set forth in the appended claims.

1. A computer implemented method for tracking a deformable global shapeof an object in motion in a sequence of digital medical images, themethod performed by the computer comprising the steps of: providing aprobabilistic subspace shape model for the deformable global shapederived from a set of training shapes; defining one or more controlpoints along the deformable global shape; tracking each of the one ormore control points as the object is in motion; estimating uncertaintyof a location of a control point in motion; exploiting the uncertaintyto constrain the deformable global shape using a point dependentcovariance matrix to employ the subspace shape constraint model using anon-orthogonal projection; fusing the subspace constraint model and theuncertainty to determine a current location of the object by correctingoutlier control points according to the probabilistic subspace shapemodel; and remedying error accumulation during tracking of the one ormore control points by using a multi-template adaptive matchingframework.
 2. The method of claim 1 wherein the step of defining the oneor more control points is performed manually.
 3. The method of claim 1wherein the step of defining the one or more control points is performedautomatically.
 4. The method of claim 1 wherein the step of trackingeach of the one or more control points employs a Bayesian kernelmatching approach.
 5. The method of claim 1 wherein the step of trackingeach of the one or more control points employs an optical flow basedapproach.
 6. The method of 1 wherein the subspace shape constraint modelis a Gaussian model.
 7. The method of claim 1 wherein the subspace shapeconstraint model is an eigenspace model.
 8. The method of claim 1wherein the object is a myocardial wall of a left ventricle.
 9. Themethod of claim 8 wherein the myocardial wall is imaged via anechocardiogram.
 10. The method of claim 8 wherein the shape tracking isused to detect abnormalities in the myocardial wall.
 11. The method ofclaim 8 wherein the myocardial wall is an endocardial wall.
 12. Themethod of claim 8 wherein the myocardial wall is an epicardial wall. 13.The method of claim 1 wherein each control point is trackedindependently of the other control points.
 14. The method of claim 1wherein the step of representing uncertainty includes measuringanisotropic uncertainty.
 15. The method of claim 1 wherein a stronglyadapted Principal Component Analysis (SA-PCA) model is used to constrainthe shape.
 16. The method of claim 15 wherein the initial shape is usedto construct the SA-PCA model.
 17. The method of claim 1 wherein theglobal shape is defined by an ordered set of control points.
 18. Themethod of claim 1 wherein the prior shape model is derived from trainingdata.
 19. The method of claim 1 wherein the shape of the object inmotion is deformable.
 20. The method of claim 19 wherein the motion isnon-rigid.
 21. The method of claim 1 further comprising the step of:displaying each subsequent shape.
 22. A program storage device readableby a computer, tangibly embodying a program of instructions executableby the computer to perform the method steps for tracking a global shapeof an object in motion in a sequence of digital medical images, methodcomprising: providing a probabilistic subspace shape model for thedeformable global shape derived from a set training shapes; defining oneor more control points along an initial contour of the deformable globalshape; tracking each of the one or more control points as the object isin motion; estimating uncertainty of a location of a control point inmotion; exploiting the uncertainty to constrain the deformable globalshape using a point dependent covariance matrix to employ the subspaceshape constraint model using a non-orthogonal projection; fusing thesubspace constraint model and the uncertainty to determine a currentlocation of the object by correcting outlier control points according tothe probabilistic subspace shape model; and remedying error accumulationduring tracking of the one or more control points by using amulti-template adaptive matching framework.
 23. The computer readableprogram storage device of claim 22 wherein the control points aredefined manually.
 24. The computer readable program storage device ofclaim 22 wherein the control points are defined automatically.
 25. Thecomputer readable program storage device of claim 22 wherein trackingeach of the one or more control points employs a Bayesian kernelmatching approach.
 26. The computer readable program storage device ofclaim 22 wherein tracking each of the one or more control points employsan optical flow based approach.
 27. The computer readable programstorage device of claim 22 wherein the subspace constraint model is aGaussian model.
 28. The computer readable program storage device ofclaim 22 wherein the subspace shape constraint model is an eigenspacemodel.
 29. The computer readable program storage device of claim 22wherein the object is a myocardial wall of a left ventricle.
 30. Thecomputer readable program storage device of claim 29 wherein themyocardial wall is imaged via an echocardiogram.
 31. The computerreadable program storage device of claim 29 wherein the shape trackingis used to detect abnormalities in the myocardial wall.
 32. The computerreadable program storage device of claim 29 wherein the myocardial wallis an endocardial wall.
 33. The computer readable program storage deviceof claim 29 wherein the myocardial wall is an epicardial wall.
 34. Thecomputer readable program storage device of claim 22 wherein eachcontrol point is tracked independently of the other control points. 35.The computer readable program storage device of claim 22 wherein theuncertainty representation includes anisotropic uncertainty.
 36. Thecomputer readable program storage device of claim 22 wherein a stronglyadapted Principal Component Analysis (SA-PCA) model is used to constrainthe shape of the contour.
 37. The computer readable program storagedevice of claim 36 wherein the initial shape is used to construct theSA-PCA model.
 38. The computer readable program storage device of claim22 wherein the global shape is defined by an ordered set of controlpoints.
 39. The computer readable program storage device of claim 22wherein the prior shape model is derived from training data.
 40. Thecomputer readable program storage device of claim 22 wherein the shapeof the object in motion is deformable.
 41. The computer readable programstorage device of claim 40 wherein the motion is non-rigid.
 42. Thecomputer readable program storage device of claim 22, the method furthercomprising the step of: displaying each subsequent shape.
 43. A computerimplemented method for tracking a deformable global shape of an objectin motion in a sequence of digital medical images, the method performedby the computer comprising the steps of: providing a probabilisticsubspace shape model for the deformable global shape derived from a setof training shapes; defining one or more control points along thedeformable global shape; tracking each of the one or more control pointsas the object is in motion; estimating uncertainty of a location of acontrol point in motion, wherein estimating uncertainty includesmeasuring anisotropic uncertainty; building multiple appearance modelsfor each control point, the multiple appearance models being based onuncertainty of a location of one or more of the control points; usingrobust information fusion to compute an estimate of the location anduncertainty of each control point; and combining the appearance models,the uncertainty, and the probabilistic subspace shape model to determinea current location of the object wherein the tracked deformable globalshape resides inside the shape model subspace and is a most plausibleshape according to uncertainties in global motion estimation local pointmeasurements, and the deformable shape model.
 44. The method of claim 43wherein the robust information fusion is Variable-BandwidthDensity-based Fusion (VBDF).
 45. The method of claim 44 wherein robustinformation fusion comprises the steps of: inputting a number ofmeasurements represented by mean vectors and covariance matrices; fusingtogether the measurements by computing variable bandwidth densityestimates; and detecting a most significant mode of the variablebandwidth density estimate.
 46. The method of claim 43 wherein theobject is a myocardial wall of a left ventricle.
 47. The method of claim46 wherein the myocardial wall is imaged via an echocardiogram.
 48. Themethod of claim 46 wherein the shape tracking is used to detectabnormalities in the myocardial wall.
 49. The method of claim 46 whereinthe myocardial wall is an endocardial wall.
 50. The method of claim 46wherein the myocardial wall is an epicardial wall.
 51. The method ofclaim 43 wherein each control point is tracked independently of theother control points.
 52. The method of claim 43 wherein a stronglyadapted Principal Component Analysis (SA-PCA) model is used to constrainthe shape.
 53. The method of claim 52 wherein the initial shape is usedto construct the SA-PCA model.
 54. The method of claim 43 wherein theglobal shape is defined by an ordered set of control points.
 55. Themethod of claim 43 wherein the shape of the object in motion isdeformable.
 56. The method of claim 55 wherein the motion is non-rigid.